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We present extensive numerical studies to determine the phase diagrams of cubic and hexagonal 
blue phases in an electric field. We confirm the earlier prediction that hexagonal phases, both 2 and 
3 dimensional, are stabilized by a field, but we significantly refine the phase boundaries, which were 
previously estimated by means of a semi-analytical approximation. In particular, our simulations 
show that the blue phase I - blue phase II transition at fixed chirality is largely unaffected by electric 
field, as observed experimentally. 
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I. INTRODUCTION 

The blue phases (BPs) of chiral molecules provide 
spectacular examples of soft solids, formed via a spon- 
taneously occurring network of disclination lines within 
a cholesteric background [l]. All BPs arise close to 
the transition between the cholesteric and the isotropic 
phase; in at least two cases (BPI and BPII) the disclina- 
tion network forms an ordered (cubic) lattice. (A third, 
BPIII, is thought to be amorphous [l].) When an elec- 
tric field is applied, several new blue phase disclination 
lattices have either been observed or theoretically pre- 
dicted. Two of these have hexagonal symmetries, and 
have been named and respectively: these have 
first been predicted theoretically [2, 3 and subsequently 
seen experimentally [4 . A further BP appearing in a 
field, BPX, has tetragonal symmetry: it was found ex- 
perimentally in the 1980s and its structure recently 
possibly revealed by simulations [6]. For a review of ex- 
perimental results on blue phases, especially in electric 
fields, the interested reader should consult Refs. [3 [8]. 

In the last few years, BPs have emerged as very promis- 
ing device materials with fast (sub-/is) and tunable color 
response. This increase in technological potential comes 
from recent remarkable experiments [Q'-^l , which man- 
aged to stabilize BPs over a temperature range of 50 K, 
compared to ranges of about 1 K previously [1 . However, 
for their potential to be fully realized, our understanding 
of BPs needs to become as quantitative as the one we 
have for conventional nematic liquid crystals. 

In particular, our theoretical knowledge of the phase 
diagrams of BP-forming liquid crystals in the presence of 
an applied field is currently quite dated. It relies on the 
seminal papers by Hornreich et al. ^ , which are how- 
ever based on an approximation (to allow semi-analytical 
progress): the tensorial liquid crystalline order parameter 
is represented by a Fourier series, comprising harmonics 
of relatively high order but only corresponding to the 
m = 2 helicity mode, which eliminates the dependence 
on one of the elastic constants (see Refs. p3 ^ for de- 
tails). These papers have provided a very useful first 
theory of blue phases, but it has become clear recently 



that their quantitative predictive power is somewhat lim- 
ited [61 [in [Ml [15] . For instance, within these approxi- 
mations, the phase diagram at zero electric field predicts 
that BPII and BPI appear in the wrong order upon vary- 
ing the molecular chirality [i2\ (cf. experimental phase 
diagrams shown in, e.g., Ref. [8]). At the same time, 
these theories were also unable to account for anoma- 
lous electrostriction of BPI [6l [161 E] • Recent simula- 
tions have shown that these shortfalls are a drawback of 
the approximations employed, and not of the underlying 
Landau - de Gennes mean field free energy [18 . The 
latter is the starting point for describing the thermody- 
namics of BPs, and when handled directly gives good 
semi-quantitative agreement with the experiments. The 
simulations of Refs. [151 [19] use a lattice Boltzmann (LB) 
method which can address the full dynamics of the sys- 
tem including fluid flow. Recent work favors instead a 
hybrid method [20l [21] in which an LB code for fluid 
degrees of freedom is coupled to a conventional finite dif- 
ference code for the order parameter (this method was 
recently used in binary fluids as well [22 ). In the current 
work, we address only thermodynamic steady states, in 
which case it is convenient to switch off the fluid motion 
and use only diffusive relaxation to find local minima of 
the free energy. (We can then compare these minima 
to construct the phase diagram.) As well as being di- 
rectly comparable with experiments, the resulting phase 
diagrams form a secure foundation for future hydrody- 
namic simulations on the same systems using the hybrid 
LB approach \2T . 

Our program in this paper is therefore to update the 
phase diagrams for blue phases in an electric field by 
computer simulation of the governing equations derived 
from the Landau - de Gennes free energy functional (de- 
tailed below). We find that the existing semi-analytical 
approximations [2l [TS] are able to capture rather 
well the qualitative physics of the problem. In partic- 
ular, we confirm that an intermediate electric field sta- 
bilizes the two previously proposed hexagonal phases, 
}i2D and HsD' However, our simulations show that the 
blue phases, whether cubic or hexagonal, extend to re- 
gions of significantly lower chirality than predicted by the 
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semi- analytical approximation; this in agreement with 
the trend which was observed in phase diagrams of BPs 
without any electric field [HI [15]. Furthermore, we also 
resolve which of the cubic blue phases, BPI or BPII, is 
stable at any given point in the phase diagram. Hence we 
find that the Landau - de Gennes free energy leads to a 
near-horizontal phase boundary in the field-temperature 
plane, in agreement with experimental observations [24l. 

It is important to note that, while in this work we focus 
on a continuum description of blue phase via a Landau - 
de Gennes free energy, there also exist simulation studies 
of individually resolved molecules, described as sphero- 
cylinders interacting with a coarse grained potential [25 , 
which managed to stabilize blue phases. This approach is 
extremely interesting, however the length and time scales 
accessible within it are significantly smaller than the ones 
we can cover with our approach, and it would be unfea- 
sible to compute a phase diagram with this method. 

This paper is organized as follows. In Section II we 
discuss the Landau - de Gennes theory which leads to 
the continuum equation of motion which we numerically 
solve. In Section III we present the numerical results, 
i.e. the phase diagrams in the chirality-temperature and 
temperature-field, planes, and we quantitatively compare 
them with the ones predicted via semi-analytical approx- 
imations. Finally, Section IV contains our discussion and 
conclusions. 



II. LANDAU - DE GENNES THEORY 
A. Equations of motion and free energy 

The nematodynamic description of liquid crystals, of- 
ten named after its inventors as the Erickson-Leslie- 
Parodi approach, uses a fixed-magnitude (headless) unit 
vector field as order parameter. This "director field" 
represents the average orientation of the liquid crystal 
molecules. However, this approach proves to be inappro- 
priate for BPs, as it cannot account for disclination lines: 
on these topological defects, the strength of nematic or- 
dering goes to zero and no average local orientation can 
be defined. A suitable description of BPs is therefore only 
possible within the framework of the Landau - de Gennes 
theory. This employs an order parameter Q, which is a 
traceless and symmetric second rank tensor; this can de- 
scribe simultaneously the direction and magnitude of lo- 
cal nematic ordering. The Q-tensor approach also allows 
for this ordering to become biaxial; in general the order 
parameter can be written as: 

Q = qi{l (g) /) + qm{rn <S) m) - ^{qi ^ g^)!, (1) 

where (g) denotes the tensorial product. The vector fields 
/ and m are two independent director fields perpendicu- 
lar to each other. The quantities qi and q^ are called the 
scalar order parameters and I is the unit tensor. This 
representation confirms Q to be a traceless symmetric 



second rank tensor which has five independent compo- 
nents. In many systems the smaller of the two scalar 
order parameters is everywhere very small (if not strictly 
zero). A director field n and a scalar order parameter q 
can then be reintroduced within the uniaxial approxima- 
tion: 

Q^q{n0n)-h. (2) 

This now has three independent components to describe 
both the strength and the director of the uniaxial order- 
ing. However, we have no need for the uniaxial approxi- 
mation in the numerical approach used in this paper. 

The phenomenological Landau - de Gennes free energy 
functional in an external electric field E is given as: 

^[Q] = |^v|^(i-|)q^^-^q«^Q;37Q7« 

+ y [iSaj6d^Qs^ + 2goQa/3)^ + {daQa^f] | (3) 

The thermodynamic equilibrium state is the global min- 
imum of for any given parameters; met ast able phases 
are local minima of J^. The complete free energy is made 
up of three different contributions. The first one (com- 
prising all terms containing Aq) is the bulk free energy, 
which contains terms of the Q-tensor up to fourth order. 
The scale factor Aq is called the bulk free energy con- 
stant, while the parameter 7 plays the role of an effective 
reciprocal temperature. For a nematogen without chi- 
rality (g'o = 0), for 7 < 2.7 the isotropic state gives the 
global minimum of the free energy, whereas for 7 > 2.7 
the system has nematic order in equilibrium. The relative 
local distortion enters the free energy functional on the 
level of first order gradient terms; throughout this paper 
we consider only the "one elastic constant" approxima- 
tion [18 , which is a common approach when investigating 
generic liquid crystallizing behavior. K is the resulting 
single elastic constant; without this approximation, Eq.[3| 
becomes considerably more complicated. We have cho- 
sen a specific representation of the gradient free energy 
yy, which incorporates a gradient-independent part of 
the bulk free energy to ensure that it is always positive. 
The parameter qo = ^ determines the intrinsic preferred 
pitch length pq of the underlying cholesteric {i.e., chiral 
nematic) liquid crystal. Thirdly, the coupling to an exter- 
nal electric field is provided by the remaining term 
in Eq. [3j with > being the dielectric anisotropy. 
This term is linear in the Q-tensor and quadratic in the 
electric field (the latter is dictated by symmetry for non- 
ferroelectric liquid crystals, as considered here). 

For the purposes of finding local minima of it is suf- 
ficient to ignore fiuid fiow and momentum conservation, 
instead taking the tensor order parameter Q to obey the 
following (purely relaxational) equation of motion: 
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The molecular field H is defined as the functional deriva- 
tive of the Landau - de Gennes free energy functional 
Eq. |3]with respect to the order parameter, and therefore 
it vanishes in equilibrium. It is specifically given by: 



^ = -;5Q + 3^''U 



(5) 



The parameter F is a rotational diffusion constant. (Note 
that in practice Eq. [4] is supplemented by an additional 
dynamical update rule for a so-called 'redshift' factor, 
described below.) 

It is convenient to render the free energy functional Eq. 
[Sjdimensionless. This gives rise to the following minimal 
set of parameters on which phase behavior can depend: 



r = 



27(1-7/3) 




327rAo7 



(6) 

(7) 
(8) 



Note that the cholesteric pitch parameter qo is used in the 
non-dimensionalization of lengths; clearly this works only 
for nonzero qo which we from now on assume. (There are 
no BPs for achiral molecules with qo = 0.) The quadratic 
term of the dimensionless bulk free energy becomes then 
proportional to the reduced temperature r, whereas the 
magnitude of the gradient free energy term is propor- 
tional to K,^ the chirality. Note that the chirality ex- 
presses a ratio between the gradient free energy and the 
bulk free energy terms and is a measure of how much 
twist the system wants to have. The parameter e is an 
effective field strength. 

The free energy functional Eq. [3] yields a highly com- 
plicated free energy landscape. Although the equilibrium 
state for each r, a^, e is uniquely defined by the functional, 
it is not generally possible to determine these states by 
any analytical minimization procedure. An exception is 
the limit of infinite chirality, where the bulk as well as the 
external field free energy terms become negligible com- 
pared to the gradient free energy and an analytical solu- 
tion can be found. While this limit of infinite chirality is 
hardly relevant for a real- world understanding or experi- 
ments, the topological character of the equilibrium defect 
structure in this limit is often unaltered upon reducing 
chirality - so long as the unit of length (related to qo) is 
rescaled appropriately. Therefore, these analytical solu- 
tions offer a useful starting point for numerical solutions 
of the relaxational dynamics: one can initialize the runs 
at finite chirality with the solution that is already known 
for the infinite chirality limit. The chosen dynamics will 
relax this structure towards a local free energy minimum; 
we then compare the minimized free energy with those 
found from competing starting structures to find the best 
possible local minimum among the ones corresponding to 
the topologies which we have considered. 



B. Initial conditions 

In this work we consider four such structures, alongside 
the isotropic phase, the (para-)nematic phase, and the 
standard uniaxial configuration for the cholesteric state: 
BPI, BPII, }i2D and Hsd (in this last case we consider 
two possible starting conditions, labeled and Hg^). 
In each case the relevant infinite-chirality expression was 
taken as initial conditions. For BPI we used: 



Qxy 
Qxz 

Qyy 

Qyz 



-2cos((7o?/) sm((7o^) + 

sin (q^x) cos (q^z) + cos (q^x) sin (q^y) 

cos {q^y) cos {q^z) + 
>/2 sin (^0^) sin (q^z) — sin (q^x) cos (q^y) 

cos((7qx) cos((7q?/) + 

sin(^o^) sm{qQy) - cos{qQx) sm{qQz) 
—2 sin(gQx) cos{qQz) + 
siniq^y) cos{qQx) + cosiq^y) siniq^z) 
V2 cos{qQz) cos{qQx) + 

V2 siniq^y) siniq^x) - siniq^y) cos(^o^), (9) 



where = V2qo. For BPII we used: 

- cos{2qoz) - cos{2qoy)] Q^y sm{2qoz) 
Qxz - sm{2qoy); Qyy cos{2qox) - cos{2qoz) 
Qyz sm{2qox). (10) 

For the hexagonal two-dimensional BP, abbreviated as 
i^2D^ we used: 

Qxx - --cos{qox)cos{V3qoy) 



Qxz 

Qyy 

Qyz 



-\/3 

— — sm{qox) sm{V3qoy) 
V3cos{qox) sm{V3qoy) 
— cos{2qox) 



cos 



{qox) cos{V3qoy) 



- sin(2g'o^) — sm{qox) cos{V3qoy). (11) 



Finally, there are two possible hexagonal three- 
dimensional BPs identified in Ref. [13] , which were found 
to be stable in the case of positive dielectric constant 
which we consider here. For the first one, abbreviated as 
Hg^, we used: 

3 1 

Qxx - --cos{qox) cos{V3qoy) + - cos{qoz) 

Vs 1 

Qxy - -^sm{qox) sm{V3qoy) + - sin(^o^) 

Qxz - cos{qox) sm{V3qoy) 

1 ^ 

Qyy - - cos{2qox) - - cos{qox) cos{V3qoy) 

Qyz - s'm{2qox) - s'm{qox) cos{V3qoy) . 



- cos{qoz) 
(12) 
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The starting condition for the second three-dimensional 
hexagonal phase, abbreviated as Hg^, is the same as the 
one used for Hg^, with the difference that all terms de- 
pending on either x or y (or both) have the opposite sign. 

In order to calculate the phase diagrams in the three- 
dimensional parameter space spanned by r, hz and , we 
follow the evolution of the free energy for the different 
configurations towards their (in general metastable) equi- 
libria by solving Eq. |4] Once the Q tensor reached a 
stationary state, we compare the free energy densities 
of the various structures and identify the lowest one as 
giving the true equilibrium. As already mentioned, this 
ignores the possibility of still lower values being achieved 
by phases that we have not considered, of which the 
most relevant are probably BPX (tetragonal) and BPIII 
(mentioned further below). In interpreting the phase di- 
agrams shown below, the reader should bear in mind 
that their validity is thus restricted to the eighth phases 
(isotropic, nematic, cholesteric, BPI, BPII, ll2D^ and the 
two phases) detailed above. We have restricted to 
these phases in line with the previous literature [2$ |13j : 
introduction of other topologies may therefore lead to 
changes in this phase diagram. 

An additional limitation to our work is as follows. We 
have already stated that the pitch parameter tz (via qo) 
sets the preferred length scale for BP structures. How- 
ever, the dimensionless number that relates the BP unit 
cell to qo is not known a priori but must be found as 
part of the free energy minimization. In numerical work, 
however, one can only simulate an integer number of unit 
cells within a periodic simulation domain. The procedure 
is then to fix the periodic domain but rescale the size of 
the physical volume that it represents by a so-called 'red- 
shift' factor. One then finds the redshift that gives the 
lowest free energy as part of the dynamical minimization. 
For a completely general disclination lattice one requires 
an independent rescaling of each spatial dimension. For 
simplicity we here assume a single redshift factor to gov- 
ern all three dimensions {i.e., we assume a 'cubic red- 
shift'). This is in line with previous literature practice 
[11 [131 [23l [2S], and it is the best we can currently do: an 
anisotropic redshift would effectively require a separate 
simulation for each aspect ratio chosen, and this would 
not be feasible with our current computational power. 
The assumption of a single redshift is exact for BPI and 
BPII only so long as these remain strictly cubic struc- 
tures. In the presence of a field, both should acquire a 
tetragonal distortion which we cannot allow for; however, 
this effect is thought to be small [2l. The two Hsd struc- 
tures are handled exactly only if the chosen aspect ratio 
of the unit cell box (which, for a cubic finite difference 
grid as used here, must be rational) coincides with the 
true one. Again, the error with our chosen value (26/15, 
see below) is expected to be minor. 

It is useful to briefly sketch how the algorithm to de- 
termine the redshift works. The method we use was orig- 
inally proposed in Ref. [15]. The free energy functional 
Eq. [3] comprises terms up to quadratic order in gradients 



of the Q-tensor: upon a rescaling of the unit cell dimen- 
sion, L L/r, it is easy to see that the gradient term 
are rescaled by a factor of r per derivative. Therefore if 
we rescale L L/r, the new free energy becomes: 



/(Q) = r^A{dQ^) + rB{dQ) + C. 



(13) 



The scaling factor r* which minimizes the functional is 
then simply 



B 

2A' 



(14) 



This observation provides a simple recipe to determine 
the redshift at every time step. 



C. Numerical Aspects 

Here we briefly present the computational details of 
our work. The size of the simulation boxes was = 



32 lattice sites for the cubic BPs BPI and 



BPII and for the cholesteric phase. To accommodate the 
hexagonal BPs we chose boxes of L^^ = 52, = = 30 
lattice sites for H2D (although a cubic box would do in 
this case) and = = 52, Ly = 30 lattice sites for the 
two HsD structures. These values were chosen as the clos- 
est approximation to the ratio 1 : >/3 which admits a 'per- 
fect' hexagonal lattice such as the HCP structure of hard 
spheres. (As mentioned above, computational limitations 
did not allow us to explore different values of this aspect 
ratio parameter, which might be somewhat different for 
the true structure, even at e = 0.) The helical pitch 
was set to 16 lattice sites for the cubic BPs and for the 
cholesteric phase, and to 15 lattice sites for the hexago- 
nal BPs. For BPI and BPII we generally simulated a box 
containing eight unit cells, to check for any (large scale) 
field-induced reconstruction; for the hexagonal BPs we 
used a single unit cell. As initial configurations we used 
a simple uniaxial helix for the cholesteric phase, together 
with the infinite chirality solutions Eqs. 9p2 The equa- 
tion of motion of the order parameter Eq.[4]was solved by 
using a finite difference scheme [27 and a rotational dif- 
fusion constant F = 0.3 in simulation units. At each time 
step, the value of the redshift, equivalent of the optimal 
scaling factor r*, was computed. This requires little extra 
calculation as the quantities A and B (see Eqs. 13-14) are 
needed for the free energy. Once we know the value of r* , 
the unit cell needs to be "rescaled" . Instead of changing 
the actual simulation cell to be simulated, which would 
be inefficient and inaccurate, we rescale the elastic con- 
stants and qo by 1/r* and r* respectively. Typical runs 
to reach equilibration required 25,000 time steps. At the 
end of the run we compared the free energy densities in 
order to determine the equilibrium phase as the one with 
the lowest free energy density. To determine the onset of 
the nematic phase a visual check of the director field was 
performed. 
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III. RESULTS 
A. Phase Diagrams 

In this Section we present the phase diagrams obtained 
for cholesteric blue phases in the presence of an electric 
field. To validate our simulations, we first checked the 
field-free case, which was published in Ref. using a 
code similar to ours (see also Ref. [14 for a version with- 
out the variable redshift). 

The left panel curves in Fig. T] were found to agree 
well with previous LB work [14] [15] except for minor 
deviations at the BPI-BPII boundary which could be ex- 
plained by a slight difference in accuracy used. A few 
general remarks are in order to relate this phase dia- 
grams to experimental ones, as well as to previous an- 
alytical predictions. At odds with the phase diagrams 
seen experimentally, such as the ones shown in [8], an- 
other BP, known as BPIII or the "Blue Fog" is missing. 
In experiments, the BPII region is entirely enclosed be- 
tween BPI and BPIII, while in Fig. [T] it emerges as a 
region which is open towards higher chiralities. BPIII 
is believed to be an amorphous BP, which is stable at 
higher chiralities and is thermodynamically distinct from 
the others [8 , 28 , 29 . Although its structure is still a sub- 
ject of discussion and has not been fully understood yet, 
there is some evidence that it might be closer to BPII 
than to the isotropic phase [8 . For theoretical attempts 
to understand the structure of BPIII see e.g. Refs. [30]- 
[32] , In principle our numerical approach does allow us 
to study large systems in regions in parameter space for 
which BPIII is expected, so that this limitation can be 
removed, at the cost of much more expensive simulations. 
(These must be many times larger than any shown here, 
so as to contain many cells of the aperiodic structure 
without strong finite size effects.) We accordingly defer 
our investigation of BPIII to a future publication [33 . 

It is interesting to compare our approach to an older 
realization of the Landau - de Gennes theory [2] (3] [12] . 
Within this approach the existence of field-stabilized 
hexagonal blue phases was first predicted. As mentioned 
in the introduction, the major drawback within this older 
framework was that, due to computational limitations at 
that time, the Q-tensor had to be expanded in a Fourier- 
series comprising several harmonics but only the m = 2 
helicity mode. This simplifications leads to the result 
that the last elastic term in Eq. (|3| does not contribute 
to the phase diagram, which we believe is at the root of 
the quantitative discrepancy with our simulations. As a 
consequence the semi-analytical theory found the incor- 
rect sequence of blue phases in the phase diagram. These 
earlier results are displayed in the right panel of Fig. [l] 
Remarkably, our approach shows instead the correct se- 
quence (cholesteric / BPI / BPII) beneath the boundary 
to the isotropic state. Furthermore, in the earlier approx- 
imate approach BPII appears adjacent to the cholesteric 
phase and BPI occurs at considerably lower temperatures 
than in the full numerical phase diagram reported here 



and in Refs. [14] [15]. Another important aspect is that 
the O5 structure, which is predicted theoretically (see 
e.g. Ref. [12]) but not observed experimentally so far, is 
relegated to higher chiralities than the ones considered 
in our phase diagrams. For a long time these issues cast 
some doubt on the appropriateness of the Landau - de 
Gennes theory for describing blue phases, but it is now 
clear that this functional is actually in good qualitative 
agreement with experiments, and that the inconsisten- 
cies with observations were essentially a drawback of the 
approximations previously employed. 

We turn now to the role of applied fields, e 7^ 0, which 
forms the main focus of the current work. Figures [2][3] 
show phase diagrams in chirality-temperature plane for 
increasing field strengths e = 0.2 and e = 0.3. We can 
identify three main effects of the external field on the 
phase boundaries. First of all, in the region where the 
isotropic phase appeared in the field-free case, the system 
becomes oriented under the influence of the external field, 
forming a paranematic phase (N) . This applies as well for 
the low-chirality part of the ordered region, where a ne- 
matic phase now appears in place of the cholesteric phase. 
For increasing field strengths the low-chirality nematic- 
cholesteric phase boundary and the cholesteric phase it- 
self both move towards the right hand side of the phase 
diagram, ie., to higher chiralities. Secondly, increasing 
the value of the electric field leaves the BPI-BPII phase 
boundary almost unchanged, an interesting feature which 
has also been observed in experiments |24j. The extent 
of the cholesteric phase is more or less retained. The BPI 
region gets therefore overtaken by the cholesteric phase 
and its phase boundaries move towards lower tempera- 
tures. Finally, very close to the BP-nematic boundary 
a pocket of hexagonal blue phases opens up. The right 
panel in Fig. [2] display where Hornreich and Shtrikman 
[2 predicted these stable hexagonal BPs at the same field 
strength. Our results again suggest that their theoreti- 
cal approach works well in terms of the general location 
of these phases, but quantitatively it suffers from the 
limitations of the approximations used. Our more accu- 
rate numerics show that while a transition from the Hsd 
phases to I12D can arise under increasing field strength 
(see Fig. [4] below), the series of transitions most often 
observed is ^sd ~^ ^20- We further predict that 

the boundary between and BPII, the closest among 
the cubic BPs, should be straight both in the (z^, r) plane. 
Finally, we confirmed the observation from earlier work 
that II2D is closer to the nematic phase with respect to 
the other hexagonal phases: this is perhaps expected as 
I12D has a symmetry plane perpendicular to the field di- 
rection and its character is therefore closer to that of a 
nematic than that of the other BPs. As the stability re- 
gion of the O5 structure is outside the chirality range we 
considered at zero field, we have not pursued its study 
further at non-zero field (we expect an electric field to 
destabilise it, as happens with the other cubic phases). 

Another way to visualize the phase behavior, espe- 
cially to determine the experimentally relevant critical 
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field strength for the switching into the nematic state, 
is to perform cuts through the r^hz^e parameter space 
for selected chiralities k. Fig. [4] and [5] show phase dia- 
grams in field-temperature plane for = 1.0 and = 1.5 
respectively. A comparison with the semi-analytical cal- 
culations of Ref. [12 shows that their phase boundaries 
are roughly in the same place as ours, except that a BPI 
phase is entirely missing in the relevant parameter region. 
The near horizontal BPI- BPII boundary seen experimen- 
tally in [24| is very well reproduced in our findings. 



B. Visualization of the BPs 

One of the most common ways to characterize the 
structure of liquid crystals is to visualize their local order, 
e.g., via their director field. Within the Q-tensor theory a 
director can be defined according to Eq. [2] as the normal- 
ized length vector which belongs to the largest eigenvalue 
of the Q-tensor. However, as the preferred configura- 
tion in BPs is an arrangement of double twist cylinders 
involving a highly complex network of disclinations, it 
is almost impossible to visualize their structure in more 
than two dimensions in this way. Perhaps the clearest al- 
ternative is then to show the defect structure itself, which 
has become the usual way to portray BP structures. In 
our simulations, this can easily be achieved by setting 
a threshold for the scalar order parameter q (defined as 
the maximum eigenvalue of the tensor order parameter) 
and imaging the resultant isosurface, which divides 'de- 
fect' (low q) from 'non-defect' (high q) zones. Although 
rather simple, this method works remarkably well and 
we use it here due to its simplicity of implementation. 
The following figures [Bp] show therefore isosurfaces of 
the scalar order parameter q. 

Fig. [6] represents a unit cell of the hexagonal two- 
dimensional BP }i2D in the field-free case e = at 
r = 0.75, = 1.5. The isosurface correspond to a value 



q = 0.11. This is a metastable configuration, because 
the II2D phase is only found to be the equilibrium phase 
for fields e ~ 0.3. Clearly recognizable is the hexagonal 
arrangement of the disclination lines. These are oriented 
along the z-direction, which is the direction of the applied 
external field. (With a field on, neither of the hexagonal 
BPs could be maintained with disclinations perpendicu- 
lar to the field.) Note that a unit cell of the hexagonal 
structures carries a full turn of the director when one 
passes through it along the x- or ^/-direction, in contrast 
with a half turn in the cases of the cubic blue phases 
(this can be seen, e.g., from checking the periodicity of 
the tensorial order parameter in Eqs. 9-12). 

If the external field is switched on, as shown in Fig. 
[7| the regions with low order become slightly larger and 
the disclination tubes facet. The magnitude of the order 
parameter at its maxima is found to be larger than in the 
field-free case. 

The situation is more complicated for the hexagonal 
three-dimensional BPs, and Hg^. Fig. [s] shows a 



metastable example of at r = 0.6, = 1.5 and e = 
0; one can readily recognize a hexagonal arrangement if 
viewing along the z direction, i.e. perpendicular to the xy 
plane. However, when viewed from other directions the 
structure looks quite intricate. The general appearance, 
but not the topology, of the structure slightly depends on 
the value of the threshold of the scalar order parameter 
q; we have chosen a threshold of q = 0.15. 

Fig. [9] shows the same phase with external field. The 
defect regions along the field direction merge and increase 
in their extent. Similarly, in the xy symmetry plane the 
domains with higher order around the center of the voids 
increase in size. 

Fig. 10 and 11 instead show the other 3D hexagonal 
phase, H^, both without and with a field. Again the for- 
mer configuration is metastable, while the latter is stable. 
With respect to Hg^, the disclination lines making up the 
Hg^ network do not merge at any point (at least for the 
values of the threshold, chosen, again q = 0.15). These 
lines are twisted and arranged onto parallel planes, the 
orientation of neighboring planes being about 60^ con- 
sistently with the overall hexagonal symmetry. When a 
field is applied and the Hg^ phase is stable, the discli- 
nations twist up (this is more evident when viewed in a 
direction perpendicular to the xy plane - compare the 
top panels of Fig. IToland 11). 



The result that the BPI-BPII phase boundary in Fig. 
|4] and [5] remains almost constant with increasing field 
strength (before the system enters first the cholesteric 
and then the nematic phase) has its counterpart in the 
near constancy of the disclination network of these struc- 
tures as the field is varied. Indeed the cubic BPs appear 
to be much more stable than the hexagonal BPs and are 
hardly deformed by the infiuence of the external field. 



The top picture in Fig. [12] gives the isosurface q = 0.22 
of the order parameter in BPI, whereas the picture at 
the bottom shows the isosurface q = 0.15 in BPII. The 
cyan surfaces depict the field-free case, while the magenta 
parts show results for the field strength e = 0.3. It is 
clear that the disclination network changes only slightly, 
and the isosurfaces are only slightly displaced in direction 
of the external field. (Recall however that our single- 
redshift approach excludes any anisotropic distortion of 
these structures that could arise at high field strengths.) 



IV. CONCLUSIONS 

We have presented extensive numerical simulations to 
investigate the phase diagram of cholesteric blue phases 
in an electric field. We have considered the case of a posi- 
tive dielectric anisotropy, and we have for computational 
simplicity taken into account only isotropic deformations 
of the unit cell for each structure. While this leaves out 
by necessity phenomena like the "anomalous electrostric- 
tion" of BPI [6l[T6l[T7], this is a good starting point, par- 
ticularly for comparisons with the earlier literature which 
often employed this approximation. 
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We confirm in tiieir main qualitative aspects the pre- 
dictions of Refs. [H [3], which were found by a semi- 
analytical approximation which led to the phase diagram 
being independent of one of the elastic constants. These 
works predicted that intermediate values of the electric 
field stabilizes two new hexagonal blue phases, one with 
2-dimensional and another one with 3-dimensional sym- 
metry. With respect to previous literature, we were able 
to refine the phase diagrams in two important ways. 
First, we find that the blue phase regions extends to sig- 
nificantly lower values of chirality than previously esti- 
mated. Second, we have established which of the cubic 
blue phases, BPI or BPII, is stable for a fixed value of 
electric field, chirality and temperature. At a given chi- 
rality value, we found that the boundary between BPI 
and BPII does not appreciably depend on temperature, 
in agreement with experimental results. Finally we have 
seen that the disclination lattices of the hexagonal blue 
phases modify quite strongly as the field is increased, 
in contrast to what happens for the cubic blue phases 
which, up until close to the transition to the nematic or 
cholesteric state, are not greatly affected in structure by 
applied fields. 

To improve further the accuracy of our phase diagrams, 
while remaining within the general framework of the Lan- 
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dau - de Gennes approach, it would be necessary not 
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This would allow further insight into the adequacy of the 
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ing for the experimental phase behavior. Meanwhile, our 
work lays a useful foundation for future computer sim- 
ulations of blue phases in electric fields, addressing for 
instance the switching dynamics of blue phase devices. 
To study dynamics, rather than the equilibrium phase 
behavior addressed here, one must however use the full 
momentum-conserving equations of motion. To this end, 
our finite difference code has already been married to 
a lattice Boltzmann approach [20l |2T] creating a hybrid 
code with which we plan to address such issues in the 
near future. 
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FIG. 1: Phase diagrams in chirality-temperature plain for e = 0, as found within our approach (left) and in Ref. [2^ (right). 
Typically, redshift values for BPI were smaller than for BPII, e.g. at r = 0, for — 0.55 and — 1.1 BPI was found to be 
stable, with a redshift r equal to, respectively, r = 0.7 and r = 0.76; while for n — lA BPII was stable and r = 0.88. 
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FIG. 2: Phase diagrams in chirality-temperature plain for e = 0.2, as found with our approach (left) and in Ref. 2 (right). 
Redshift values were typically smaller for BPI than for BPII, and smaller for cubic blue phases than for hexagonal ones: e.g. 
for BPI at r = 0, and hi = 1.1, r = 0.82; for BPII at r = and At = 1.4, r = 0.89; for Hsz? at r = 0.95 and hi = 0.9, r = 0.96; 
for at r = 0.9 and Hi = 0.9, r = 0.935; for at t = 0.8 and Hi — 0.9, r = 0.91. Note that for smaller values of e we 
reproduce electrostriction and the redshift decreases albeit very slightly, in agreement with previous theoretical and numerical 
literature p[23]. 
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FIG. 3: Phase diagram in chirality-temperature plain for e = 0.3. 
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FIG. 4: Phase diagrams in field strength-temperature plain for k = 1.0 as found with our approach (left) and in Ref. [2 (right). 
Note that change in scale for r. 
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FIG. 5: Phase diagrams in field strength-temperature plain for k = 1.5 as found with our approach (left) and in Ref. [2^ (right). 




FIG. 6: Hexagonal-2D BP in equilibrium state at r 
0.75, HI — 1.5, e = 0. Each tube contains a disclination line. 




FIG. 7: Hexagonal-2D BP in equilibrium state at r = 
0.75, = 1.5, e = 0.3. 




FIG. 8: Hexagonal-3Da BP in metastable state at r 
0.6, = 1.5, e = 0. 




FIG. 9: Hexagonal-3Da BP in equilibrium state at r = 
0.5, = 1.25, e = 0.3. 




FIG. 10: Hexagonal-3Db BP in metastable state at r = 
0.5, K= 1.25, e = 0. 




FIG. 11: Hexagonal-3Db BP in equilibrium state at r = 
0.6, AC = 1.5,6 = 0.3. 




FIG. 12: BPI (top, at T = -0.8, = 1.5) and BPII (bottom, 
at T = 0, At = 1.5) at e = (cyan) and e — 0.3 (magenta). 



